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Abstract. This paper develops the theory of abehan Routh reduction for discrete mechanical 
systems and applies it to the variational integration of mechanical systems with abelian 
symmetry. The reduction of variational Runge-Kutta discretizations is considered, as well as 
the extent to which symmetry reduction and discretization commute. These reduced methods 
allow the direct simulation of dynamical features such as relative equilibria and relative 
periodic orbits that can be obscured or difficult to identify in the unreduced dynamics. 

The methods are demonstrated for the dynamics of an Earth orbiting satellite with a 
non-spheiical J2 correction, as well as the double spherical pendulum. The J2 problem is 
interesting because in the unreduced picture, geometric phases inherent in the model and those 
due to numerical discretization can be hard to distinguish, but this issue does not appear in 
the reduced algorithm, where one can directly observe interesting dynamical structures in the 
reduced phase space (the cotangent bundle of shape space), in which the geometric phases 
have been removed. 

The main feature of the double spherical pendulum example is that it has a nontrivial 
magnetic term in its reduced symplectic form. Our method is still efficient as it can directly 
handle the essential non-canonical nature of the symplectic structure. In contrast, a traditional 
symplectic method for canonical systems could require repeated coordinate changes if one is 
evoking Darboux' theorem to transform the symplectic structure into canonical form, thereby 
incurring additional computational cost. Our method allows one to design reduced symplectic 
integrators in a natural way, despite the noncanonical nature of the symplectic structure. 



1. Introduction 

This paper addresses reduction theory for discrete mechanical systems with abehan symmetry 
groups and its relation to variational integration. To establish the setting of the problem, a few 
aspects of the continuous theory are first recalled (see |29 1 for general background). 

Continuous Reduction Theory. Consider a mechanical system with configuration manifold 
Q and a symmetry group G (with Lie algebra g) acting freely and properly on Q and 
hence, by cotangent lift on T*Q, with corresponding (standard, equivariant) momentum 
map J : T*Q — > g*. Recall from reduction theory (see ^36' "25'], and references therein) 
that, under appropriate regularity and nonsingularity conditions, the flow of a G-invariant 
Hamiltonian H : T*Q ^ M. naturally induces a Hamiltonian flow on the reduced space 
= 3^^{iJ.)/Gf^, where is the isotropy subgroup of a chosen point /i G g*. In the 
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abelian case, if one chooses a connection 21 on the principal bundle Q Q/G, then is 
symplectically isomorphic to T*{Q/G) carrying the canonical symplectic structure modified 
by magnetic terms; that is terms induced from the /i-component of the curvature of 21. 

The Lagrangian version of this theory is also well-developed. In the abelian case, it goes 
by the name of Routh reduction, (see, for instance |29|, §8.9). The reduced equations are 
again equations on T{Q /G) and are obtained by dropping the variational principle, expressed 
in terms of the Routhian, from Q to Q/G. The nonabelian version of this theory was originally 
developed in I32II33I . with important contributions and improvements given in I15ll30i . 

Of course, reduction has been enormously important for many topics in mechanics, such 
as stability and bifurcation of relative equilibria, integrable systems, etc. We need not review 
the importance of this process here as it is extensively documented in the Uterature. 

Purpose, Main Results, and Examples. This paper presents the theory and illustrative 
numerical implementation for the reduction of discrete mechanical systems with abelian 
symmetry groups. The discrete reduced space has a similar structure as in the continuous 
theory, but the curvature will be taken in a discrete sense. The paper studies two examples 
in detail, namely satellite dynamics in the presence of the bulge of the Earth (the J2 effect) 
and the double spherical pendulum (which has a nontrivial magnetic term). In each case the 
benefit of studying the numerics of the reduced problem is shown. Roughly, the reduced 
computations reveal dynamical structures that are hard to pick out in the unreduced dynamics 
in a way that is reminiscent of the phenomena of pattern evocation, as in 1341 1351 . Another 
interesting application of the theory is that of orbiting multibody systems, studied in |42 41 1. 

We refer to 1 37 1 for a review of discrete mechanics, its numerical implementation, some 
history, as well as references to the literature. The value of geometric integrators has been 
documented in a number of references, such as |9|. In the present paper, we shall focus, to 
be specific, on discrete Euler-Lagrange and variational symplectic Runge-Kutta schemes and 
their reductions. One could, of course, use other schemes as well, such as Newmark, Stormer- 
Verlet, or Shake schemes. However, we wish to emphasize that without theoretical guidelines, 
coding algorithms for the reduced dynamics need not be a routine procedure since the reduced 
equations are not in canonical fonn because of nontrivial magnetic terms. For example, using 
Darboux' theorem to put the structure into canonical form so that standard algorithms can be 
used is not practical. We also remind the reader that there are real advantages to taking the 
variational approach to the construction of symplectic integrators. For example, as in 1231 . 
the variational approach provides the design flexibility to take different time steps at spatially 
different points in an asynchronous way and still retain all the advantages of symplecticity 
even though the algorithms are not strictly symplectic in the naive sense; such an approach is 
well-known to be useful in molecular systems, for instance. 

Motivation for Discrete Reduction. Besides its considerable theoretical interest, there are 
several practical reasons for carrying out discrete Routh reduction. These are as follows: 

(i) Features that are clear in the reduced dynamics, such as relative equilibria and 
relative periodic orbits, can be obscured in the unreduced dynamics, and ap- 
pear more complicated through the process of reconstruction and associated ge- 
ometric phases. This is related to the phenomenon of pattern evocation that 
is an important practical feature of many examples, such as the double spher- 
ical pendulum 1341 135 1 and the stepping pendulum |12|. Going to a suitable 
(but non-obvious) rotating frame can "evoke" such phenomena (see the movie at 

|http : //www . cds . caltech . edu7~marsden7research/demos/movies/Wendlandt /pattern . mpq) . 
This is essentially a window to the reduced dynamics, which the theory in the present 
paper allows one to compute directly. 
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(ii) While directly studying the reduced dynamics can yield some benefits, it can be difficult 
to code using traditional methods. In particular, the presence of magnetic terms in the 
reduced symplectic form, as is the case with the double spherical pendulum, means 
that traditional symplectic methods for canonical systems do not directly apply; if one 
attempts to do so, it may result (and has in the literature) in many inefficient coordinate 
changes when evoking Darboux' theorem to put things into canonical form. 

(iii) Although simulating the reduced dynamics involves an initial investment of time in 
computing geometric quantities symbolically, these additional terms do not appreciably 
affect the sparsity of the system of equations to be solved. As such, direct coding of the 
reduced algorithms can be quite efficient, due to its reduced dimensionality. 

Two Obvious Generalizations. The free and proper assumption that we make on the group 
action means that we are dealing with nonsingular, that is, regular reduction (see [381 for the 
general theory of singular reduction and references to the literature). It would be interesting to 
extend the work here to the case of singular reduction but already the regular case is nontrivial 
and interesting. While our examples have singular points and the dynamics near these points 
is interesting, there is no attempt to study this aspect in the present paper. 

Secondly, it would be interesting to generalize the present work to the case of nonabelian 
groups and to develop a discrete version of nonabelian Routh reduction, (as in 1 15 30|). We 
believe that such a generalization will require the further development of the theory of discrete 
connections, which is currently part of the research effort on discrete differential geometry 
(see |,20|, and references therein.) Other future directions are discussed in the conclusions. 

Other Discrete Reduction Results. We briefly summarize some related results that have been 
obtained in the area of reduction for discrete mechanics. First of all, there is the important 
case of discrete Euler-Poincare and Lie-Poisson reduction that were obtained in l2l 1271 l28l . 
This theory is appropriate for rigid body mechanics, for instance. 

Another important case is that of discrete semidirect product reduction that was obtained 
in (3] El and applied to the case of the heavy top, with interesting links to discrete elastica. 
This case is of interest in the present study since the heavy top, as with the general theory of 
semidirect product reduction (see II31II13I ). one can view the reduction of this problem as 
Routh reduction. Linking these two approaches is an interesting topic for future research. 

Outline. After recalling the notation from continuous reduction theory, fj2]develops discrete 
reduction theory, derives a reduced variational principle, and proves the symplecticity of 
the reduced flow. The relationship between continuous- and discrete-time reduction is also 
discussed. How the variational (and hence symplectic) Runge-Kutta algorithm induces a 
reduced algorithm in a natural way is shown in fj3l In 50] we put together in a coherent 
way the main theoretical results of the paper up to that point. In ^the numerical example 
of satellite dynamics about an oblate Earth is given, and in fj6l the example of the double 
spherical pendulum, which has a non-trivial magnetic term, is given. Lastly, in we address 
some computational and efficiency issues. 

2. Discrete Reduction 

In this section, it is assumed that the reader is familiar with continuous reduction theory as 
well as the theory of discrete mechanics; reference is made to the relevant parts of the literature 
as needed. It will be useful to first recall some facts about discrete mechanical systems with 
symmetry (see, for instance, |37 1 for proofs.) 
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Discrete Mechanical Systems with Symmetry. Let G be a Lie group (which shortly will be 
assumed to be abelian) that acts freely and properly (on the left) on a configuration manifold 
Q. Given a discrete Lagrangian : Q y. Q ^ M. that is invariant under the diagonal action 
of G on Q X Q, the corresponding rfiscreie momentum map : Q x Q ^ g* is defined by 

^d{qQ,qi) ■ ? = D2Ld{qo,qi) ■ (.giqi), (i) 

where D2 denotes the derivative in the second slot and where is the infinitesmal generator 
associated with ^ G g. The map is equivariant with respect to the diagonal action of G on 
Q X Q and the coadjoint action on g*. The discrete Noether theorem states that the discrete 
momentum is conserved along solutions of the DEL (Discrete Euler-Lagrange) equations, 

D2Ldiqk-i,qk) + DiLd{qk,qk+i) =0. (2) 

Notice that 

Jd('7o,9i) • C = J (£"2 id (90,91)) • 

where J : T*Q g* is the momentum map on T*Q; i.e., Jd = J o Fi^, where 
VLd ~ D2Ld : Q X Q ^ T*Q is the discrete Legendre transform. Thus, for fi e g*, 
we have FL^ (Jj^^(/x)) C J^^(/i). The symplectic algorithm (usually called the position- 
momentum form of the algorithm) obtained on T*Q from that on Q x Q via the discrete 
Legendre transform thus preserves the standard momentum map J. 

There will be a standing assumption in this paper, namely that the given discrete 
Lagrangian Ld is regular; that is, for a point {q,q) € Q x Q on the diagonal, the iterated 
deiiwatiwe D2DiLd{q, q) : TqQxTqQ R is a non-degenerate bilinear form. By the implicit 
function theorem, this implies that a point {qk-i, qk) near the diagonal and the DEL equations 
(|3 uniquely determines the subsequent point qk+i in a neighborhood of the diagonal in Q x Q 
(or, if one prefers, for small time steps); in other words, the DEL algorithm is well defined 
by the DEL equations. Regularity also implies that the discrete Legendre transformation 
¥Ld = D2Ld : Q X Q T*Q is a local diffeomorphism from a neighborhood of the 
diagonal in Q x Q to a neighborhood in T*Q. For a detailed discussion, see |37 1. 

Reconstruction. The following lemma gives a basic result on the reconstruction of discrete 
curves in the configuration manifold Q from those in shape space, defined to be S* = Q/G. 
The Lemma is similar to its continuous counterpart, as in, for example, 1151 . Lemma 2.2. The 
natural projection to the quotient will be denoted ttq.g : Q Q/G; q ^ x = [q\G (the 
equivalence class of q g Q). Let YeT{q) denote the vertical space at q; namely the space of 
all vectors at the point q that are infinitesimal generators S,Q{q) G TqQ; or in other words, the 
tangent space to the group orbit through q. We say that the discrete Lagrangian Ld is group- 
regular if the bilinear map D2DiLd{q, q) : TqQ x TqQ M restricted to the subspace 
Yer{q) x YeT{q) is nondegenerate. In addition to regularity, we shall make group-regularity 
a standing assumption in the paper The following result is fundamental for what follows. 

Lemma 2.1 (Reconstruction Lemma). Fix /i e g* and let xo,xi, ... ,Xn be a sufficiently 
closely spaced discrete curve in S. Let qo, qi £ Q be such that [go]G = a^o, [qi\G — a^i cind 
Jdiqo, 9i) = A*- Then there is a unique closely spaced discrete curve qi, q2, . . . , qn such that 
[qk]G = Xk and Jd{qk-i,qk) = l^Jor k = l,2,...,n. 

Proof. We must construct a point q2 close to qi such that [q2]G = X2 and Jd{qi, 92) = /i; the 
construction of the subsequent points (73, . . . g„ then proceeds in a similar fashion. 

To do this, pick a local trivialization of the bundle ttq g : Q ^ Q/G where locally 
Q — S X G and write points in this trivialization as qo = (xq, go) and qi = {xi,gi), etc. 
Given the points go — (2^0,50), qi — {xi,gi) with Jd(go, qi) — we seek a near identity 
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group element k ^ G such that 172 ■— {x2,kgi) satisfies Jd(<Zi, 92) = According to 
equation Q, this means that we must satisfy the condition -D2^ci('?ij 92) • Cq(92) = (Mj ; 
for all ^ e g. In the local trivialization, this reads 

D2Ld{{xi,gi), {x2,kgi)) ■ {Q,TRkg,(.) = , 0) 

where Rg is right translation on G by 5. Consider solving the equation 

D2Ldi{xi,gi),ix2,kgi)) ■ iO,TRkg,0 = {^^,0 , (4) 

for fc as a function of the variables gi,xi, X2 with /i fixed. By assumption, there is a solution 
for the case xi ~ xo, X2 — xi and gi ~ ga, namely k ~ ko ^ 9i9o^ "ear identity group 
element). The implicit function theorem shows that when the point go,xo,xi is replaced by 
the nearby point gi,xi, X2, there will be a unique solution for k near fco provided that the 
derivative of the defining relation (|3j with respect to k at the identity is invertible, which is 
true by group-regularity. Since group regularity is G-invariant, the above argument remains 
valid as fc; drifts from the identity. □ 

Note that the above Lemma makes no hypotheses about the sequences Xn or q„ satisfying 
any discrete evolution equations. 

To carry out reconstruction in the continuous case, in addition to the requirements that 
the lifted curve in TQ lie on the momentum surface, and that it projects to the reduced 
curve x{t) e S under ttq g, one also requires that it be second-order, which is to say that 
it is of the form {q{t),q{t)). If a connection is given, then the lifted curve is obtained by 
integrating the reconstruction equation-again, see flS] for details. The discrete analogue 
of the second-order curve condition is explained as follows. Consider a given discrete 
curve as a sequence of points, (0:0, xi), {xi, X2), • ■ • , {xn-i, Xn) in S* x 5. Lift each of the 
points in S* X S* to the momentum surface J^^(/i) C Q x Q. This yields the sequence, 
(^Q, g"), {ql,ql), . . . , ((70^^,(7"^^), which is unique up to an overall diagonal group action. 
The discrete analogue of the second-order curve condition is that this sequence in Q x Q 
defines a discrete curve in Q, which corresponds to requiring that qi — q^^^, for k ~ 
0, . . . , n — 1, which is clearly possible in the context of the reconstruction lemma. 

Discrete reconstruction naturally leads to issues of discrete geometric phases, and it 
would be interesting to express the discrete geometric phase in terms of the discrete curvature 
on shape space; this will surely involve some ideas from the currently evolving subject of 
discrete differential geometry and so we do not attempt to push this idea further at this point. 

While many of the computations we present in this paper are in the setting of local 
trivializations, the results are valid globally through the construction given below. 

Identification of tlie Quotient Space. Now assume that G is abelian so that G — 
acts on J^^(//) and so that the quotient space J^^(/x)/G makes sense. We assume the 
above regularity hypotheses and freeness and properness of the action of G so that this 
quotient is a smooth manifold. It is clear that the map ipf^ : J'^^{iJ.)/G —> S x S given 
by [{'l,q')]G {x,x') is well-defined, where the square brackets denotes the equivalence 
class with respect to the given G action, and where x = [qlc- The argument given in the 
reconstruction lemma shows that for a point (qa, qi) & J^^ (/i), (/?^ is a local diffeomorphism 
in a neighborhood of the point [(qoi 9i)]g- In fact, the uniqueness part of that argument shows 
that for two nearby points (qi, 92) and {q[, q!^) in J^^(/x), if qi = giq[ and q2 = .9292' '^^en 
91 =92- Thus, there is a neighborhood [/ of a given a chain of closely spaced points lying in 
Jj^dj,) with this property. Saturating this neighborhood with the group action, we can assume 
that U is G-invariant. Restricted to U, ip^ becomes a diffeomorphism to a neighborhood of 
the diagonal of 5 x 5*. 
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Assume, as above, that L^; : Q x Q ^ M is a discrete Lagrangian that is invariant under 
the action of an abehan Lie group G on Q x Q. In view of the preceding discussion, L^i 
restricted to J^^ (/i) (and in the neighborhood of a given chain of points in this set) induces a 
well defined function Ld{xQ,xi) of pairs of points (a;o, xi) m S x S. This discrete reduced 
Lagrangian will play an important role in what follows. 

Discrete Reduction. Let q := {go, • • • , In] be a solution of the discrete Euler-Lagrange 
(DEL) equations. Let the value of the discrete momentum along this trajectory be ^. Let 
Xi — [<li\G^ SO that X := {xq, . . . , Xn] is a discrete shape space trajectory. Since q satisfies the 
discrete variational principle, it is appropriate to ask if there is a reduced variational principle 
satisfied by x. 

An important issue in dropping the discrete variational principle to the shape space 
is whether we require that the varied curves are constrained to he on the level set of 
the momentum map. The constrained approach is adopted in [151, and the unconstrained 
approach is used in L30.I . In the rest of this section, we will adopt the unconstrained approach 
of 1301 . and will show that the variations in the discrete action sum, without assuming the 
variations at the endpoints vanish and evaluated on a solution of the discrete Euler-Lagrange 
equations, depends only on the quotient variations, and therefore drops to the shape space 
without constraints on the variations. 

If q is a solution of the DEL equations, the interior terms in the variation of the discrete 
action sum vanish, leaving only the boundary terms; that is, 

ri-l 

6 ^ Ld{qk, qk+i) = DiLd{qo, qi) ■ Sqo + D2Ld{qn-i,qn) ■ Sqn- (5) 

fc=0 

Given a principal connection 21 on Q, there is a horizontal-vertical spUt of each tangent 
space to Q denoted Vq = horvq + veiVq for Vq e TqQ. Thus, 

D2Ld{qn-i,q7i) ■ Sqn = D2Ld{qn-i,qn) ■ hor(5q„ + D2Ld{q„-i,q„) ■ veiSqn- 

As in continuous Routh reduction, we will rewrite the terms involving vertical variations 
using the fact that we are on a level set of J^;. Namely, write the vertical variation as 
verSqn — £,Q{qn), where ^ ~ 2t(5(7„) and use the definition Q of to give 

D2Ld(qn-i,q7i) ■ verSqn = D2Ld{qn-i, qn) ■ Cqiqu) = Jd(gn-i,gn) • C 

(o) 

= = (M,21((5(?„)) = %{qn) ■ Sqn- 

Thus, the boundary terms can be expressed as 

D2Ld{qn-i-,qn) ■ 5qn = D2Ld{qn-i, qn) ' hor(5q„ + 2lp(q„) • Sqn, (7) 
DiLdiqo,qi) ■ 6qo = DiLd{qo,qi) -horSqo -'Qif^iqo) ■ Sqo, (8) 

and so the variation of the discrete action sum, becomes 

ri-l 

S ^ Ld{qk, qk+i) = DiLd{qo, qi) ■ horSqo + D2Ld{qn-i,qn) ■ hor5g„ 

+ ■ Sqn - 2l,,(<Zo) • Sqo, (9) 

when restricted to solutions of the discrete Euler-Lagrange equations. 

Motivated by the preceding equation, introduce the 1-form ^ on Q x Q defined by 

A = 7r*2%, - 7rl%, (10) 
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where 7ri,7r2 : Q x Q Q are projections onto the first and the second components 
respectively. This allows us to expand the boundary terms involving 21^ into a telescoping 
sum, and rewrite (|9} in terms of the 1-form A as 

ji-i 

^{DLd - A)(qk,qk+i) ■ i^Qk^Sqk+i) = DiLa{qo,qi) ■ horSqo 

+ D2Ld{qn-l,qn) -^OTSqn- (11) 

We now drop il 1\ to the reduced space SxS. Consider the projection maps tt : QxQ ^ {Qx 
Q)/G, and tt^^^ : ^(m) ~^ ^(/^)/^' ^"'^ inclusion maps, ; J^^(/i) ^ Q x Q, 
and ln^d '■ Jrf ^ (m) /G ^ [Q x Q)/G. Then clearly, the following diagram commutes, 

3-\fi)/G(^^iQxQ)/G 

By the G-invariance, Ld drops to a function Ld on the quotient {Q x Q)/G, so that 
Ld — LdO TT. The pullback (in this case, the restriction) of Ld to J^^(/i)/G' is called the 
discrete reduced Lagrangian and is denoted Ld- Thus, Ld = LdoZ^^d', identifying ^ (/i) /G 
with SxS, this definition agrees with Ld as defined earlier 

Lemma 2.2. T/ie 1-fonn AonQxQ restricted to 3'^^ (/i), drops to a 1-form A on ^ (/i) /G 
anc/ induces, (for closely spaced points), via the map Lp^^, a 1-form that we denote by the same 
letter, on S x S. Similarly, the 1-form DLd ^ Aon Q x Q restricted to the momentum level 
set J^^ (/i) then drops to the 1-form DLd — Aon on J^^ (/x) / G and induces, via the map (/?^, 
a 1-form that we denote by the same letter, on SxS. 

Proof. The equivariance of the projections tt^ with respect to the diagonal action on Q x Q 
and the given action on Q, together with the invariance of the 1-form 21^ on Q imply that A 
is invariant. Since 2t^ vanishes on vertical vectors for the bundle Q Q/G, it follows that 
A vanishes on vertical vectors for the bundle Q x Q ^ [Q x Q)/G. Therefore, there is a 
1-form ^ on (Q X Q) /G such that A = n*A. 

Since Ld = LdO tt, and the exterior derivative commutes with pull-back, it follows that 
dLd = n*dLd. From tt o ^ = l^^d o Tr^.d, we get t* = = '^l.Jl.A-^- Thus, 

the 1-form A restricted to J^^(/i) drops to the 1-form A = Z*^ ^.4 on 3'^^{p)/G. Similarly, 
i* ^did = 71"* ^t* ddLd and so dLd restricted to J^^(/i) drops to the 1-form t* ^jdZ^ = dLd 
on3~^^{p)/G. These 1 -forms push forward under the map (^^ : 3'^^{fj,)/G —^Sx 5* in the 
manner that was explained earlier □ 

With the preceding lemma, and equation il 1\ . we conclude that 

^{DLd - A){xk,Xk+i) ■ {Sxk,Sxk+i) = DiLd{qo,qi) ■ \\or5qo 

k—0 

+ D2Ldiqn-i,qn) ■ horSqn- (12) 

Assuming that Sx vanishes at the endpoints, hor Sqo = 0, and hor 5qi = and consequently, 
the boundary terms vanish and we obtain the reduced discrete variational principle, 

n—1 n—1 

6^Ld{xk,Xk+i) = ^A(xk,Xk+i) ■ {Sxk,6xk+i). (13) 

fe=0 fc=0 
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In an analogous fashion to rewriting DLii{xk, Xk+i)-{Sxic, Sx^+i) as DiLd{xk,Xk+i)-6xk + 
D2Lii{xk, Xk+i) ■ Sxk+i, we do the same by for the A term by defining, 

A{xo,xi) ■ {Sxo,Sxi) = ^i(a;o,a;i) • Sxo + ^2(2:0, a^i) • Sxi. 

Then, equating terms involving Sxk on the left hand side of ( I13> to the corresponding terms 
on the right, we get the discrete Routh (DR) equations giving dynamics on 5 x S*: 

D2Ld{xk-i,Xk)+DiLdixk,Xk+i) = A(xfe_i,a;fc)+yli(xfe,a;fc+i).(14) 

Note that these equations depend on the value of momentum /it. Thus, if q is a discrete curve 
satisfying the discrete Euler-Lagrange equations, the curve x obtained by projecting q down 
to S satisfies the DR equations (I14> . 

Now consider the converse, the discrete reconstruction procedure: Given a discrete curve 
X on S* that satisfies the DR equations, is x the projection of a discrete curve q on Q that 
satisfies the DEL equations? 

Let the pair (go, 9i) be a lift of {xo,xi) such that Jd(<7o, 9i) ~ fJ- Let q = {go, • ■ • , Qn} 
be the solution of the DEL equations with initial condition (go, 91)- Note that q has 
momentum /i. Let x' = {x'q, . . . , xJj} be the curve on S obtained by projecting q. By our 
arguments above, x' solves the DR equations. However x' has the initial condition {xq, xi), 
which is the same as the initial condition of x. By uniqueness of the solutions of the DR 
equations, x' = x. Thus x is the projection of a solution q of the DEL equations with 
momentum /i. Also, for a given initial condition qq, there is a unique lift of x to a curve 
with momentum fi. Such a lift can be constructed using the method described in Lemma IZTI 
Thus, lifting x to a curve with momentum /i yields a solution of the discrete Euler-Lagrange 
equations, which projects down to x. 

We summarize the results of this section in the following Theorem. 

Theorem 2.3. Let x is a discrete curve on S, and let q be a discrete curve on Q with 
momentum fi that is obtained by lifting x. Then the following are equivalent. 

1 . q solves the DEL equations. 

2. q is a solution of the discrete Hamilton's variational principle 

ri-l 
k=0 

for all variations (5q o/q that vanish at the endpoints. 

3. X solves the DR equations 

D2Ld{xk^i,Xk) + DiLd{xk,Xk+i) = A2{xk-i,Xk) + Ai{xk, Xk+i). 

4. X is a solution of the reduced variational principle 

n—1 n—1 

d^Ldixk,Xk+i) = ^A{xk,Xk+i) ■ {Sxk,Sxk+i) 

k=a k=0 

for all variations Sx ofx that vanish at the endpoints. 

Note that for smooth group actions, the order of accuracy will be equal for the reduced 
and unreduced algorithm. 



Discrete Routh Reduction 



9 



2.1. Preservation of the Reduced Discrete Symplectic Form 

The DR equations define a discrete flow map Fk : S x S ^ S x S. We already know that 
the flow of the DEL equations preserves the symplectic form onQ x Q. In this section 
we show that the reduced flow F^ preserves a reduced symplectic form fi^ on S* x S", and 
that this reduced symplectic form is obtained by restricting fi^^ to ^ (/i) and then dropping 
to 5* X 5*. In other words, tt* j^^f_i,d — i*^ rf^id- The continuous analogue of this equation is 

Recall from continuous reduction theory on the Hamiltonian side that in the case of 
cotangent bundles, the projection tt^ : J^^(/i) — > T*S can be defined as follows: If 
aq € J~^{ii), then the momentum shift aq — 2lp(g) annihilates all vertical tangent vectors 
at (J e Q, as shown by the following calculation: 

{aq - %{q),^Q{q)) = J{aq) . e - (a^, = (m, - (A^, = 0. 
Thus, aq — 2lp((7) induces an element of T*S and 7r^(Q!g) is defined to be this element. 



Let F' : J^^(m) ^ J^^ip) be the restriction of FL^ to J;^^(Ai). Thus F' 



o t 



d o ¥Ld, where : J^^{ii) T*Q and ^ : J^^(Ai) — > Q x Q are inclusions. Define 
the map F : S* x 5 ^ T*S by F(a;o,a;i) = D2Ld{xo,xi) - ^2(2:0, a;i). By equation 
and Lemma lz!2l this map is well-defined. The map F will play the role of a reduced discrete 
Legendre transform, and in contrast to the continuous theory, the momentum shift appears in 
the map F, as opposed to the projection 7r^,d. As in the continuous theory, ftfj,^d does involve 
magnetic terms. 

Lemma 2.4. The following diagram commutes. 



S xS — -^T*S 

Proof Let (170,91) e J^^m)- Thus D2Ld{qo,qi) e J^^ifJ.), and 
7rp(F'((?o, 91)) = '^tj,{D2Ld{qa, qi))- 

As noted above, tt^ (£'2^^(90, 91)) is the element of T*_^S determined by 
{D2Ld{qo, qi) - %{qi))- For 6qi £ Tq^Q, we have 

{D2Ldiqo,qi)~'^M'Sll) = {DLd - A){qo,qi) ■ {0,6qi) 

= (DLd- A)ixo,xi) ■ {0,6x1) 
= D2Ld{xQ,xi) -Sxi -A2ixo,xi) ■ Sxi , 

where in the second equality, we used Lemma l2?2l Thus, 

T^p.{D2Ld{qa,qi)) = D2Ld{xo,xi) - A2{xa,xi), 
which means F o vr^^rf = tt^ o F'. □ 
Theorem 2.5. The flow of the DR equations preserves the symplectic form 

^tj.,d = F*(r25 - 7r|,.5 

where [3^ is the 2-form on S obtained by dropping dSt^. 

Furthermore, fi^.d can be obtained by dropping to S x S the restriction of Qlci to 



^(/Lt). In other words, TT^^j^^fx,d = i^*^,d^ 



Li- 
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Proof. We give an outline of the steps involved; the details are routine to fill in. The strategy 
is to first show that the restriction to Jj^(/x) of the symplectic form fi^^^ drops to a 2-form 
ri^ tj on 5 X S*. The fact that the discrete flow on Q x Q preserves the symplectic form $7^^ 
is then used to show that the reduced flow preserves 51 ^t ^. The steps involved are as follows. 

(i) Consider the l-formG^^ on Q x Q defined by 0l^((7o, 9i) • {5qa,5qi) D2Ld{qo, qi)- 
Sqi. The 1-form 6^^ is G-invariant, and thus the Lie derivative C^^^i^Ql^ is zero. 

(ii) Since fl^^ — —dQ^^, fl^^ is G-invariant. If d : J^^(/i) — > Q x Q is the inclusion, 
®Ld ~ 'ui.d^Li and 17'^^ = i*^^jP-Li are the restrictions of Ql^ and VLl^ respectively to 
J (/i). One checks that and fi'^^ are invariant under the action of G on J^^ (/it). 

(iii) If is an infinitesimal generator on Jj^^(/i), then 

^..-(.) ^ - -^..-(.) ^ del. - -^«.-u.Qi^. + ^ Q'-. - 0- 

This follows from the G-invariance of B'^^, and the fact that 8^^ • Cj-i(^) = (Mj 0- 

(iv) By steps 2 and 3, the form 51'^^ drops to a reduced form fi^.^j on J^^(/i)/G ~ S x S. 
Thus, if TT^.d : J^^(/x) — > x 5* is the projection, then tt* j^^^^d — ^'l^- Note that the 
closure of il^ ^ follows from the fact that is closed, which in turn follows from the 
closure of fi^^ and the relation il^^ = t* ^^La- 

(\) If Fk : Q X Q ^ Q X Q is the flow of the DEL equations, let be the restriction 
of this flow to J^^ (/i). We know that drops to the flow Fk of the DR equations on 
S X S. Since Fk preserves il^^, F^ preserves 51/^^. Using this, it can be shown that F/. 
preserves fi^ rf. Note that it is sufficient to show that tt* ^(-F^fl^.d) — ti"^ d^t^-d- 
(vi) It now remains to compute a formula for the reduced form fi^^d. Using Lemma 1241 it 
follows that 

Thus TT* rfflp,(i = TT* ^F*(ri5 - TTy.g g/?^), from which it follows that fi^^d F*(fi5 - 
TT^.g s/^Ai)- Incidentally, this expression shows that 57^ ^ is nondegenerate provided the 
map F = id — -^2 is a local diffeomorphism. 

□ 

One can alternatively prove symplecticity of the reduced flow directly from the reduced 
variational principle — see §2.3.4 of |20|. 

2.2. Relating Discrete and Continuous Reduction 

As shown in [371, if the discrete Lagrangian Ld approximates the Jacobi solution of the 
Hamilton-Jacobi equation, then the DEL equations give us an integration scheme for the EL 
equations. In our commutative diagrams we will denote the relationship between the EL and 
DEL equations by a dashed arrow as follows: 

{TQ,EL) - - ^{Q X Q,DEL). 

This arrow can thus be read as "the corresponding discretization". By the continuous and 
discrete Noether theorems, we can restrict the flow of the EL and DEL equations to J^^{^) 
and J^^(/i) respectively. The flow on Jj^^{ij) induces a reduced flow on J£^{ii)/G ~ TS, 
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which is the flow of flie Routh equations. Similarly, the discrete flow on J^^(/i) induces 
a reduced discrete flow on J^^(/Lt)/G w 5 x S*, which is the flow of the discrete Routh 
equations. Since the DEL equations give us an integration algorithm for the EL equations, it 
follows that the DR equations give us an integration algorithm for the Routh equations. 
To numerically integrate the Routh equations, we have two options: 

(i) First solve the DEL equations to yield a discrete trajectory on Q, which can then be 
projected to a discrete trajectory on S. 

(ii) Solve the DR equations to directly obtain a discrete trajectory on Q. 

Either approach yields the same result, which we express by the commutative diagram: 



{Jl\^i),EL) 



(TS, R) ■ 



{3-^\li),DEL) 



{S X S, DR) 



(15) 



3. Reduction of the symplectic Runge-Kutta algorithm. 



A well-studied class of numerical schemes for Hamiltonian and Lagrangian systems is the 
symplectic partitioned Runge-Kutta (SPRK) algorithms (see 1 10 1 1 1 for history and details). 
We will adopt a local trivialization to express the SPRK method in which the group action 
is addition. Given a connection 2t on Q, it can be represented in local coordinates as 
21(6', x){9, x) — A{x)x + 9. By rewriting the symplectic partitioned Runge-Kutta algorithm 
in terms of this local trivialization, and using the local representation of the connection, we 
obtain the following algorithm on T*S: 



Xi = Xq + 



Si = So 



dA 



Xi = Xq + h y ^ 



Si = So + 



^R^' 
dx 

dRt" 
dx 



{Xj,Xj) 



dA 



inA{Xi) - ^lA{xo)) 



(16fl) 
(16^) 
(16c) 
(160 

(I6e) 
(16/) 



This system of equations is called the reduced symplectic partitioned Runge-Kutta (RSPRK) 
algorithm. A detailed derivation can be found in §2.5 of |20|. As we obtained this system by 
dropping the symplectic partitioned Runge-Kutta algorithm from J^^(/i) to T*S, it follows 
that this algorithm preserves the reduced symplectic form il^ ^ fig — tt^.^ on T*S. 

Since the SPRK algorithm is an integration algorithm for the Hamiltonian vector field 
Xh on T*Q, the RSPRK algorithm is an integration algorithm for the reduced Hamiltonian 
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vector field on T*S. The relationship between cotangent bundle reduction and the 
reduction of the SPRK algorithm can be represented by the following commutative diagram: 

(J-1(m),Xh) - - - {J-\ti),SPRK) 



{T*S,Xh,) 



{T*S, RSPRK) 



The dashed arrows here denote the corresponding discretization, as in jl5l l. The SPRK 
algorithm can be obtained by pushing forward the DEL equations by the discrete Legendre 
transform. See, for example, 1 37 1 . By Lemma lZ4l this implies that the RSPRK algorithm can 
be obtained by pushing forward the DR equations by the reduced discrete Legendre transform 
¥ ^ D2Ld — A2- These relationships are shown in the following commutative diagram: 



{3-/{^Ji),DEL)■ 



(5 X S, DR) 



{J-\n),SPRK) 



{T*S, RSPRK) 



4. Putting Everything Together 



The relationship between Routh reduction and cotangent bundle reduction can be represented 
by the following commutative diagram: 



{JZ\fi),EL) 



(TS, R) ■ 



¥L 



{J-^{fi),XH) 



{T*S,Xh, 



We saw in fj2j2|that if i^j approximates the Jacobi solution of the Hamilton-Jacobi equation, 
discrete and continuous Routh reduction is related by the following diagram: 

{JZ\^Ji),EL)--^{J-\^i),DEL) 



{TS,R)- 



{S X S, DR) 



The dashed arrows mean that the DEL equations are an integration algorithm for the EL 
equations, and that the DR equations are an integration algorithm for the Routh equations. 

Pushing forward the DEL equation using the discrete Legendre transform ¥Ld yields the 
SPRK algorithm on T*Q, which is an integration algorithm for X^- This is depicted by: 



{Jl\p),EL) 



iJ^\f,),DEL) 



Fid 



{J-\^,),XH)- 



- ^ {J-^ in), SPRK) 



Discrete Routh Reduction 



13 



J-\^i),XH 



J-Hp),SPRK 




Jl\pi),EL- 




-3-/{^i),DEL 



T*S,Xh^ 




■T*S,RSPRK 




TS,R X S,DR 



Figure 1. Complete commutative cube. Dashed arrows represent discretization from tlie 
continuous systems on the left face to the discrete systems on the right face. Vertical an'ows 
represent reduction from the full systems on the top face to the reduced systems on the bottom 
face. Front and back faces represent Lagrangian and Hamiltonian viewpoints, respectively. 



The SPRK algorithm on J~^{n) C T*Q induces the RSPRK algorithm on J-'^{^j.)/G « 
T*S. As we saw in ^j3] this reduction process is related to cotangent bundle reduction and to 
discrete Routh reduction as shown in the following diagram: 



iJ-Hp),XH) 



{T*S,Xh,) 



{J-\fi),SPRK) 



{T*S, RSPRK) 



iJ^Hp),DEL) 



{S X S, DR) 



Putting all the above commutative diagrams together into one diagram, we obtain Figure^ 



5. Example: J2 Satellite Dynamics 
5.7. Configuration Space and Lagrangian 

An illustrative and important example of a system with an abelian symmetry group is that of 
a single satellite in orbit about an oblate Earth. The general aspects and background for this 
problem are discussed in 1 40 1 , and some interesting aspects of the geometry underlying it, are 
discussed in [7J. 

The configuration manifold Q is R'^, and the Lagrangian is 



L{qA)^\Ms\\qf 



MsViq), 



where Ms is the mass of the satellite and : M'^ ^ M is the gravitational potential due to the 
earth truncated at the first term in the expansion in the ellipticity 



Viq) = 
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Here, G is the gravitational constant, is the mass of the Earth, is the radius of the 
earth, J2 is a small non-dimensional parameter describing the degree of ellipticity, and is 
the third component of q. In non-dimensional coordinates, 

L{qA)^\M 

This corresponds to choosing space and time coordinates in which the radius of the Earth is 1 
and the period of orbit at zero altitude is 2tt when J2 = (spherical earth). 



J2 



2 llal|2 



(17) 



5.2. Symmetry Action 

The symmetry of interest to us is that of rotation about the vertical (q^) axis, so the symmetry 
group is the unit circle S^. Using cylindrical coordinates q — (r, 6*, z) for the configuration, 
the symmetry action is : {r,9,z) 1-^ {r,9 + (f>,z). Since \\q\\, \\q\\, and g"^ = 2: are all 
invariant under this transformation, so too is the Lagrangian. 

This action is clearly not free on all of Q = IR"^^ as the z-axis is invariant for all group 
elements. This is not a serious obstacle as the lifted action is free on T{Q \ (0, 0, 0)) and this 
is enough to permit the application of the intrinsic Routh reduction theory. Alternatively, one 
can simply take Q = M^^ \ {(0, 0, z) | z G M} and then the theory literally applies. 

The shape space S = Q/G is, thus the half-plane S — x R and we will take 
coordinates (r, z) on 5*. In doing so, we are implicitly defining a global diffeomorphism 
SxG ^ Q given by {{r,z),e) ^ {r.O.z). 

The Lie algebra q for G = S*^ is the real line g = M, and we will identify the dual 
with the real line itself g* = E. For a Lie algebra element ^ G g, the corresponding 
infinitesimal generator is given by S,q : {r,9,z) ((r, 6', z), (0, ^, 0)). The Lagrange 
momentum map Jl : TQ g* is given by 3L{vq) ■ ^ = {¥L{vq),£^Q{q)), which in 
our case is a scalar quantity, the vertical component of the standard angular momentum: 
JL{ir,9,z),if,e,z)) = rH. 

Consider the Euclidean metric on R-^, which corresponds to the kinetic energy norm in 
the Lagrangian. From this metric we define the mechanical connection 21 : TQ g given 
by 2t((r, 9, z), (r, 9, z)) = 9. The 1-form 21^ on Q is thus given by 21^ = iid9. The exterior 
derivative of this expression gives d2l^ = pi(P9 = 0, and so the reduced 2-form is f3^ ~ 0. 



5.3. Equations of Motion 

The Euler-Lagrange equations for the Lagrangian ( I17> gives the equations of motion. 



J2 



3(g^) 



3\2 



To calculate the reduced equations, we begin by calculating the Routhian 

i?'^(r, 9, z, f, 9, z) = L(r, 9, z, r, 9, z) - 2l^(r, 9, z) ■ (r, 9, z) 

J2 fiz^ 



2r2 



- ^i9. 



We choose a fixed value /i of the momentum and restrict ourselves to the space ^(z^), on 
which 9 = ji. The reduced Routhian : TS ^ R is the restricted Routhian dropped to the 
tangent bundle of the shape space. In coordinates this is 



R^{r,z,f,z) = \\\{i-,z)f- 



J2 



3z2 
2r2 



1 2 
2^- 
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Recalling that (3^ — 0, the Routh equations can be evaluated to give 



"1 




/3z2 




r 




v2r2 


-01 



which describes the motion on the shape space. 

To recover the unreduced Euler-Lagrange equations from the Routh equations one uses 
the procedure of reconstruction. This is covered in detail in I26ll25ll30l . 

5.4. Discrete Lagrangian System 

We discretize this system with a high-order discrete Lagrangian. Recall that the pushforward 
discrete Lagrange map associated with this discrete Lagrangian is a symplectic partitioned 
Runge-Kutta method. 

Given a point {qo,qi) E Q x Q we will take {qo,po) and (<7i,pi) to be the associated 
discrete Legendre transforms. As the discrete momentum map is the pullback of the canonical 
momentum map, we have that JL^(go, qi) = {pe)o = {pe)i - Take a fixed momentum map 
value /i and restrict Ld to the set J£^^(/i). Dropping this to S* x S* now gives the reduced 

discrete Lagrangian : 5 x S* ^ R. More explicitly, Ld depends on {rk,0k,Zk) and 
{rk+i,Ok+i, Zk+i), but group invariance implies that the group variables only enters in the 
combination {9k+i — Ok). The condition 3d{qk, qk+i) = can be inverted to eliminate the 
dependence of Ld on {0k+i — Sk)^ and Ld can be expressed in terms of /i, (r^, zj.), and 
{rk+i,Zk+i), which yields Ld. 

As discussed in fj3] the fact that we have taken coordinates in which the group action is 
addition in 6 means that the pushforward discrete Lagrange map associated with the reduced 
discrete Lagrangian is the reduced method given by (I16a|16/t . In fact, as the mechanical 
connection has A{r, z) = and /3p ~ 0, the pushforward discrete Lagrange map is exactly a 
partitioned Runge-Kutta method with Hamiltonian equal to the reduced Routhian. These are 
generically related by a momentum shift, rather than being equal. 

Given a trajectory of the reduced discrete system, we can reconstruct the unreduced 
discrete trajectory by solving for 6. Correspondingly, a trajectory of the unreduced discrete 
system can be projected onto shape space to give a trajectory of the reduced discrete system. 

5.5. Example Trajectories 

We compute the unreduced trajectories using the fourth-order SPRK algorithm, and the 
reduced trajectories using the corresponding fourth-order RSPRK algorithm. 

Solutions of the Spherical Earth System. Consider initially the system with J2 — 0. This 
corresponds to the case of a spherical Earth, and so the equations reduce to the standard 
Kepler problem. A slightly inclined circular trajectory is shown in Figure |2 in both the 
unreduced and reduced pictures. Note that the graph of the reduced trajectory is a quadratic, 
as = vr^ + z^ is a constant. 

We will now investigate the effect of two different perturbations to the system, one due 
to taking non-zero J2 and the other due to the numerical discretization. 

The J2 effect. Taking J2 = 0.05 (which is close the actual value for the Earth), the system 
becomes near-integrable and experiences breakup of the KAM tori. This can be seen in Figure 
|3] where the same initial condition is used as in Figure|2l 

Due to the fact that the reduced trajectory is no longer a simple curve, there is a 
geometric -phase-like effect which causes precession of the orbit. This precession can be seen 
in the thickening of the unreduced trajectory. 
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Figure 2. Unreduced (left) and reduced (right) views of an inclined elliptic trajectory for the 
continuous time system with J2 = (spherical Earth). 



0.4 




Figure 3. Unreduced (left) and reduced (light) views of an incHned elHptic trajectory for the 
continuous time system with J2 = 0.05. Observe that the non-spherical terms introduce 
precession of the near-elliptic orbit in the symmetry direction. 

Solutions of the Discrete System for a Spherical Earth. We now consider the discrete 
system with J2 = 0, for the second order Gauss-Legendre discrete Lagrangian with stepsize 
ofh = 0.3. The trajectory with the same initial condition as above is given in Figure|4] 

As can be seen from the reduced trajectory, the discretization has caused a similar 
breakup of the periodic orbit as was produced by the non-zero J2. This induces precession 
of the orbit in the unreduced trajectory, in a way which is difficult to distinguish from the 
perturbation above due to non-zero J2 if only the unreduced picture is considered. If the 
reduced pictures are consulted, however, then it is immediately clear that the system is much 
closer to the continuous time system with J2 = than to the system with non-zero J2. 

Solutions of the Discrete System with J2 Effect. Finally, we consider the discrete system 
with non-zero J2 — 0.05. The resulting trajectory is shown in Figure |5l and it is clearly 
not easy to determine from the unreduced picture whether the precession is due to the J2 
perturbation, the discretization, or some combination of the two. 

Taking the reduced trajectories, however, immediately shows that this discrete time 
system is structurally much closer to the non-zero J2 system than to the original J2 = 
system. This confusion arises because both the J2 term and the discretization introduce 
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Figure 4. Unreduced (left) and reduced (right) views of an inclined trajectory of the discrete 
system with step-size h = 0.3 and J2 = 0. The initial condition is the same as that used 
in figure I2] The numerically introduced precession means that the unreduced picture looks 
similar to that of figure|3]with non-zero J2, whereas only by considering the reduced picture 
we can see the correct resemblance to the J2 = case of figure I2I 
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Figure 5. Unreduced (left) and reduced (right) views of an inclined trajectory of the discrete 
system with step-size h = 0.3 and J2 = 0.05. The initial condition is the same as that used 
in figure|3] The unreduced picture is similar to that of both figures|3]and|l] By considering the 
reduced picture the correct resemblance tol3l 

perturbations which act in the symmetry direction. 

While this system is sufficiently simple that one can run simulations with such small 
timesteps that the discretization artefacts become negligible, this is certainly not generally 
possible. This example demonstrates how knowledge of the geometry of the system can 
be important in understanding the discretization process, and how this can give insight into 
the behavior of numerical simulations. In particular, understanding how the discretization 
interacts with the symmetry action is extremely important. 

5.6. Coordinate Systems 

In this example we have chosen cylindrical coordinates, thus making the group action addition 
in 6. One can always do this, as an abelian Lie group is isomorphic to a product of copies 
of M and S^, but it may sometimes be preferable to work in coordinates in which the group 
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action is not addition. For example, cartesian coordinates in the present example. Reasons 
for choosing a different coordinate system might include ease of computation, or simplicity 
of the expressions. 

If we adopt a coordinate system wherein the group action is not expressed in terms of 
addition, the RSPRK method is not applicable, but we can still apply the Discrete Routh 
equations to obtain an integration scheme on S* x 5*. The push forward of this under F yields 
an integration scheme on T* S. The trajectories on the shape space that we obtain in this 
manner could be different from those we would get with the RSPRK method. However, in 
both cases we would have conservation of symplectic structure, momentum, and the order of 
accuracy would be the same. One could choose whichever approach is cheaper and easier. 

6. Example: Double Spherical Pendulum 

6.1. Configuration Space and Lag rang ian 

We consider the example of the double spherical pendulum which has a non-trivial magnetic 
term and constraints. The configuration manifold Q is 5^ x 5^, and the embedding linear 
space F is R'^ X K.^. The position vectors of each pendulum with respect to the pivot point are 
denoted by qi and q2. These vectors are constrained to have lengths li and I2 respectively, 
and the pendula masses are denoted by mi and 1712- The Lagrangian is 

-^-(qi,q2,qi,q2) imi||qi|p + ^TO2||qi +q2|p - migqi ■k-m2g{qi + CI2) ■ k, 

where g is the gravitational constant, and k is the unit vector in the z direction. The constraint 
function c : V ^ M.^ is given by, c(qi,q2) — (|lqi|| — h, |lq2|| — h)- Using cylindrical 
coordinates q.; = (r^, 9i, Zi), L becomes 



-migzi - m2g (21 + 22) , 

where ip — 62 — 9i. Furthermore, we can automatically satisfy the constraints by performing 
the substitutions, Zi = {If — )^/^, and Zi = —rifi{yjl'j — rf)~^/'^. 

6.2. Symmetry Action 

The symmetry of interest to us is the simultaneous rotation of the two pendula about 
vertical (z) axis, so the symmetry group is the unit circle S^. Using cylindrical coordinates 
qi = (r^, 6*^, Zj) for the configuration, the symmetry action is : {ri,di,Zi) 1-^ {ri,9i + (j), Zi). 
Since ||qi||, ||qi||, ||qi + q2||> and qi ■ k are all invariant under this transformation, so too is 
the Lagrangian. 

This action is clearly not free on all of y = x R'^, as the z-axis is invariant for all group 
elements. However, this does not pose a problem computationally, as long as the trajectories 
do not pass through the downward hanging configuration, corresponding to ri = r2 = 0. 
To treat the downward handing configuration properly, we would need to develop a discrete 
Lagrangian analogue of the continuous theory of singular reduction described in |38 1. 

We will now show the geometric structures involved in implementing the RSPRK 
algorithm. See §2.10.2 of |20| for a detailed derivation. The Lie algebra g for G = 5^ is 
the real line g = M, and we will identify the dual with the real line itself g* = R. For 
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a Lie algebra element ^ G g, the corresponding infinitesimal generator is given by : 

{ri,Bi,zi,r2,02,Z2) ^ ((ri, 6li, zi, ra, 6*3, Z2), (0, ^, 0, 0, ^, 0)). As such, the momentum 
map is given by 

JL{{ri,0i,Zi,r2,92, Z2), {riJi, Zi,r2,02, h)) 

= (toi + 7712) rjOi + m2rl92 + m2rir2 [Oi + ^2) cosip + {rif2 - r2'ri) simp, 

which is simply the vertical component of the standard angular momentum. 
The locked inertia tensor is given by 1251 . 

I(qiq2) = rnillqj^ll^ + TO2||(qi + q2)^||^ = mirf + ni2{r\ + + 2rir2 cosip). 

The mechanical connection as a 1-form is given by 

a(qi, q2) = [rriirl + ■m2{r\ + + 2rir2 cosip)]~^ [{mi + 1712) rfdOi + m2r2d02 

+m2rir2 (dOi + 6.62) cos ip + (ridr2 — r2dri) sin ip] . 

The /x-component of the mechanical connection is given by, 

a/x(qi,q2) = ^J■[mlrl + m2{rl + r| + 2rir2 cos.^)]"^ 

X { [("^1 + ^2) rf + m2rir2 cosip] d9i + [m2rl + m2rir2 coa(p] 662} . 

Taking the exterior derivative of this 1-form yields a non-trivial magnetic term on the reduced 
space, which drops to the quotient space to yield, 

= ^1712 [2 (mi + 1112) rir2 + (rnir^ + ni2{rl + r|)) cosip] 

X [mir^ + TO2 {1^1 + f 2 + 2rir2 cos 95)] ^ d(p A (r2dri — ridr2). 

The local representation of the connection is given by 



A{ri,r2, ip) — m2{mirl + m2{r\ + r| + 2rir2 cosp)) ^ 



—r2 sm(/3 
ri sin Lp 
+ rir2 cos Lp 



and the amended potential has the form 

VM = -mig{ll~rlY/^-m2g 



{I 



2)1/2 



-2 ^ii^[mirl+m2{rl+rl + 2rir2COSip)] ^ 



y^ . Recall that 



The Routhian on the momentum level set is given by, = ^ ||hor(g, v) 
hor(ug) — Vq — ^Q{vq), whcrc ^ = a{vq), and ^Q{vq) = (0, 0, 0, ^, 0). Then we obtain, 

hor(z;,) = - (0, 0, 0, a(uq), 0) = (ri,^i - ii, 7^2, ^2 - a(w,),i2). 

The kinetic energy metric has the form 



mi + m2 








m2 cos p> 


—171272 sin ip 








(mi + m2)r^ 





m2ri sin Lp 


m2rir2 cos ip 











mi + m2 











m2 cos ip 


m2ri sin ip 





m2 








m2r2 siniy9 


m2rir2 cos ip 








m2r| 




















m2 



This together with the expression for hor(uq) allows us to compute i||hor(q, and when 
combined with the formula for the amended potential gives the Routhian i?^. Notice that 
all our expressions are in terms of the reduced variables on TS, so they trivially drop to yield 
i?^. These expressions can then be directly substituted into the RSPRK algorithm in equations 
J16a|16/) to obtain the example trajectories presented in the next subsection. 
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6.3. Example Trajectories 

We have computed the reduced trajectory of the double spherical pendulum using the fourth- 
order RSPRK algorithm on the Routh equations. We first consider the evolution of ri, r2, 
and ip, using the RSPRK algorithm on the Routh equations, as well as the projection of the 
relative position of 1112 with respect to nii onto the xy plane as seen in Figure|6l 

FigureQillustrates that the energy behavior of the trajectory is very good, as is typical of 
variational integrators, and does not exhibit a spurious drift. In contrast, the non-symplectic 
fourth-order Runge-Kutta applied to the unreduced dynamics has a systematic drift in the 
energy, even when using time-steps that are smaller by a factor of four 
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Figure 6. Time evolution of ri, r2, f, and the trajectory of m2 relative to mi using RSPRK. 
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Figure 7. Relative energy drift {E — Eq)/ Eq using RSPRK (left) compared to the relative 
energy drift in a non-symplectic RK (right). 
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7. Computational Considerations 

Reduced vs. Unreduced Simulations. The reduced dynamics can either be computed 
directly, by using the discrete Routh or RSPRK equations, or by computing in the unreduced 
space, and projecting onto the shape space. We discuss the relative merits of these approaches. 

Given a configuration space and symmetry group of dimensions n and m, respectively, 
we can either use a simpler algorithm in 2n dimensions, or a more geometrically involved 
algorithm in 2(n — m) dimensions. The reduced algorithm involve curvature terms that need 
to be symbolically precomputed, but these do not affect the sparsity of the system of equations, 
and the lower dimension results in computational saving that are particularly evident in the 
repeated, or long-time simulation of problems with a shape space of high codimension. 

An example which is of current engineering interest is the dynamics of connected 
networks of systems with their own internal symmetries, such as coordinated clusters of 
satellites modeled as rigid bodies with internal rotors. If the systems to be connected are 
all identical, the geometric quantities that need to be computed, such as the mechanical 
connection, have a particularly simple repeated form, and the small additional upfront effort 
in implementing the reduced algorithm can result in substantial computational savings due to 
the lower dimension of the reduced system. 

Intrinsic vs. Non-intrinsic methods. Non-intrinsic numerical schemes, such as SPRK 
applied to the classical Routh equations, can have undesirable numerical properties due to the 
need for coordinate dependent local trivializations and the presence of coordinate singularities 
in these local trivializations, as is the case when using Euler angles for rigid body dynamics. 

In the case of non-canonical symplectic forms, frequent computationally expensive 
coordinate changes are necessary when using standard non-intrinsic schemes, as documented 
in 1441 B9I. due to the need to repeatedly apply Darboux' theorem to put the symplectic 
structure into canonical form. In contrast, intrinsic methods do not depend on a particular 
choice of coordinate system, and allow for the use of global charts through the use of 
containing vector spaces with constraints enforced using Lagrange multipliers. 

Coordinate singularities can affect the quality of the simulation in subtle ways that may 
depend on the numerical scheme. In the simulation of the double spherical pendulum, we 
notice spikes in the energy corresponding to times when ri or r2 are close to 0. These errors 
accumulate in the non-symplectic method, but remain well-behaved in the symplectic method. 
Alternatively, sharp spikes can be avoided altogether by evolving the equations on a constraint 
surface in M.^ x R-^, as opposed to choosing local coordinates that automatically satisfy the 
constraints. Here, the increased cost of working in the containing linear space with constraints 
is offset by not having to transform between charts of Sf_^ x Sf^ , which can be significant 
if the trajectories are chaotic. An extensive discussion of the issue of representations and 
parametrizations of rotation groups and its implications for computation can be found in |5 1. 

Methods which exhibit local conservation properties on each chart, may still exhibit a 
drift in the conserved quantity across coordinate changes. As discussed in [T], only methods 
in which the local representatives of the algorithm commute with coordinate changes exhibit 
geometric conservation properties that are robust. In particular, intrinsic methods retain their 
conservation properties through coordinate changes. 

Another promising approach is to use the exponential map to update the numerical 
solution, which is the basis of Lie group integrators, see 1141 . The Lie group approach 
can yield rigid body integrators that are embedded in the space of 3 x 3 matrices but 
automatically evolve on the rotation group, without the use of constraints or reprojection. 
In I16II17I . analogues of the explicit Newmark and midpoint Lie algorithm are presented that 
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automatically stay on the rotation group, and exhibit good energy and momentum stability. 
A method based on generating functions and the exponential on Riemannian manifolds is 
introduced in 1 19 1. In the variational setting, a Lie group variational integrator for rigid body 
dynamics is described in L18J . and higher-order generaUzations can be found in L21J . 

8. Conclusions and Future Work 

This paper derives the Discrete Routh equations on the discrete shape space S x S, which 
are symplectic with respect to a non-canonical symplectic form, and retains the good energy 
behavior typically associated with variational integrators. Furthermore, when the group 
action can be expressed as addition, we obtain the Reduced Symplectic Partitioned Runge- 
Kutta algorithm on T*S, that can be considered as a discrete analogue of cotangent bundle 
reduction. By providing an understanding of how the reduced and unreduced formulations are 
related at a discrete level, we enable the user to freely choose whichever formulation is most 
appropriate, and provides the most insight to the problem at hand. 

Certainly one of the obvious things to do in the future is to extend discrete reduction to 
the case of nonabelian symmetry groups following the nonabelian version of Routh reduction 
given in fl5' TO]. There are also several problems, including the averaged J2 problem, in 
which one can carry out discrete reduction by stages and relate it to the semidirect product 
work of This is motivated by the fact that the semidirect product reduction theory of fT3l 
is a special case of reduction by stages (at least without the momentum map constraint), as 
was shown in [6J . 

In further developing discrete reduction theory, the discrete theory of connections on 
principal bundles developed in [22] is particularly relevant, as it provides an intrinsic method 
of representing the reduced space {Q x Q)/G as {S x S) (B G, where G = Q Xq G with 
G acting by conjugation on G. Indeed, such a construction can be viewed as a connection on 
a Lie groupoid, and it is natural to express discrete mechanics on Q x Q in the language of 
pair groupoids, as originally proposed in [43J. Generalizations of this approach to arbitrary 
Lie groupoids, as well as a discussion of the role of discrete connections in yielding a discrete 
analogue of the Lagrange-Poincare equations, can be found in [24]. 

Another component that is needed in this work is a good discrete version of the calculus 
of differential forms. Note that in our work we found, being directed by mechanics, that the 
right discrete version of the magnetic 2-form is the difference of two connection 1-forms. 
It is expected that we could recover such a magnetic 2-form by considering the discrete 
exterior derivative of a discrete connection form in a finite discretization of space-time, and 
taking the continuum limit in the spatial discretization. Developing a discrete analogue of 
Stokes' Theorem would also provide insight into the issue of discrete geometric phases. Some 
preliminary work on a discrete theory of exterior calculus can be found in fsl . 
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